#include <deque>
#include <algorithm>

#include "ctreecontourextractor.h"
#include "augmentedctree.h"
#include "edgecache.h"

using namespace PBVP;

CTreeContourExtractor::CTreeContourExtractor(Reg3DVolumeData* pData, AugmentedContourTree* pTree)
{
	m_pData = pData;
	m_pVolCritical = new Volume3DCritical(pData);
	m_pConTree = pTree;

	m_pData->getDimension(m_Dim);
	m_pData->getSpan(m_Span);
	m_pData->getOrigin(m_Orig);

	m_Inner = 0;
}

CTreeContourExtractor::~CTreeContourExtractor()
{
	delete m_pVolCritical;
}

GeometryData* CTreeContourExtractor::extractContour(Reg3DVolumeData* pReg3D, float isoval)
{
	GeometryData *pSurf = new GeometryData();

	std::vector<SuperArc>* pArcVector = m_pConTree->cut(isoval);
	for (int i = 0; i < pArcVector->size(); i++)
	{
		int startTetra = getStartTetra((*pArcVector)[i], isoval);
		GeometryData *pContour = extractContourComponent(isoval, startTetra);
		pSurf->merge(pContour);
		delete pContour;
	}
	delete pArcVector;
	return pSurf;
}

int CTreeContourExtractor::extractLargestContours(float isoval,
												   std::vector<GeometryData*> &contours,
												   int nc /* = 1 */)
{
	std::vector<Contour3DData*> convector;
	std::vector<SuperArc>* pArcVector = m_pConTree->cut(isoval);
	for (int i = 0; i < pArcVector->size(); i++)
	{
		int startTetra = getStartTetra((*pArcVector)[i], isoval);
		Contour3DData *pContour = extractContourComponent(isoval, startTetra);
		convector.push_back(pContour);
		// printf("contour %d volume = %f, area = %f\n", i, pContour->volume(), pContour->area());
	}
	delete pArcVector;

	std::sort(convector.begin(), convector.end(), CompareContour());
	int count;
	int vecsize = convector.size();
	for (count = 0; count < nc && count < vecsize; count++)
	{
		printf("extracted contour %d: volume = %f, surface area = %f\n", 
			count, convector[vecsize-count-1]->innerVolume(), 
			convector[vecsize-count-1]->area());
		contours.push_back(convector[vecsize-count-1]);
	}

	return count;
}

Contour3DData* CTreeContourExtractor::extractContourComponent(float isoval, int startTetra)
{
	Contour3DData* pGeoData = new Contour3DData();

	// total number of tetrahedron = 6 * number of cells
	int numTetra = 6 * m_pData->numOfCells();
	bool *touched = new bool[numTetra];
	memset(touched, 0, numTetra);

	// The queue of tetrahdron to be marched
	std::deque<int> tetraQueue;
	// add the first tetra
	touched[startTetra] = true;
	tetraQueue.push_back(startTetra);

	// edge cache
	EdgeCache cache;

	float face_rang[2];

	while (!tetraQueue.empty())
	{
		int tid = tetraQueue.front();
		tetraQueue.pop_front();

		float vals[4], x0[3], x1[3], x2[3], x3[3];
		m_pVolCritical->getTetraData(tid, vals, x0, x1, x2, x3);

		Tetra *tetra = new Tetra(x0, x1, x2, x3, vals[0], vals[1], vals[2], vals[3]);

		int vids[4];
		m_pVolCritical->getTetraVertIds(tid, vids);
		// extract isosurf triangles from the tetra and add them to the surface
		tetra->extractSurface(vids, isoval, pGeoData, cache);
		pGeoData->addArea(tetra->isosurfAera(isoval));
		if(m_Inner == 0) {
			pGeoData->addInnerVolume(tetra->innerVolume(isoval));
		} else {
			pGeoData->addInnerVolume(tetra->getVolume() - tetra->innerVolume(isoval));
		}
		delete tetra;
		// consider four face-neighbor tetrahedra
		for (int i = 0; i < 4; i++)
		{
			face_rang[0] = std::min(vals[(i+1)%4], std::min(vals[(i+2)%4], vals[(i+3)%4]));
			face_rang[1] = std::max(vals[(i+1)%4], std::max(vals[(i+2)%4], vals[(i+3)%4]));

			if ((face_rang[0] <= isoval && face_rang[1] >= isoval)) // tetra intersets the contour
			{
				int neighbor_tid = m_pVolCritical->tetNeighbor(tid, i);
				if (neighbor_tid >= 0 && !touched[neighbor_tid])
				{
					touched[neighbor_tid] = true;
					tetraQueue.push_back(neighbor_tid);
				}
			}
		}
	}
	delete[] touched;
	// Compute isosurface normal vectors using function gradients
	calcSurfNormals(pGeoData);
	return pGeoData;
}

Contour3DData* CTreeContourExtractor::extractContourComponentWithVolume(float isoval, int startTetra)
{
	Contour3DData* pGeoData = new Contour3DData();

	// total number of tetrahedron = 6 * number of cells
	int numTetra = 6 * m_pData->numOfCells();
	bool *touched = new bool[numTetra];
	memset(touched, 0, numTetra);

	// The queue of tetrahdron to be marched
	std::deque<int> tetraQueue;
	// add the first tetra
	touched[startTetra] = true;
	tetraQueue.push_back(startTetra);

	// edge cache
	EdgeCache cache;

	float face_rang[2];

	while (!tetraQueue.empty())
	{
		int tid = tetraQueue.front();
		tetraQueue.pop_front();

		float vals[4], x0[3], x1[3], x2[3], x3[3];
		m_pVolCritical->getTetraData(tid, vals, x0, x1, x2, x3);

		Tetra *tetra = new Tetra(x0, x1, x2, x3, vals[0], vals[1], vals[2], vals[3]);

		int vids[4];
		m_pVolCritical->getTetraVertIds(tid, vids);
		// extract isosurf triangles from the tetra and add them to the surface
		tetra->extractSurface(vids, isoval, pGeoData, cache);
		pGeoData->addArea(tetra->isosurfAera(isoval));
		if(m_Inner == 0) {
			pGeoData->addInnerVolume(tetra->innerVolume(isoval));
		} else {
			pGeoData->addInnerVolume(tetra->getVolume() - tetra->innerVolume(isoval));
		}
		delete tetra;
		// consider four face-neighbor tetrahedra
		for (int i = 0; i < 4; i++)
		{
			face_rang[0] = std::min(vals[(i+1)%4], std::min(vals[(i+2)%4], vals[(i+3)%4]));
			face_rang[1] = std::max(vals[(i+1)%4], std::max(vals[(i+2)%4], vals[(i+3)%4]));
			
			if ((face_rang[0] <= isoval && face_rang[1] >= isoval) // tetra intersets the contour
				|| (m_Inner == 0 && face_rang[1] < isoval ) // tetra belongs to inner region of the contour
				|| (m_Inner == 1 && face_rang[0] > isoval ) )
			{
				int neighbor_tid = m_pVolCritical->tetNeighbor(tid, i);
				if (neighbor_tid >= 0 && !touched[neighbor_tid])
				{
					touched[neighbor_tid] = true;
					tetraQueue.push_back(neighbor_tid);
				}
			}
		}
	}
	delete[] touched;
	pGeoData->setOuterVolume(m_pData->getVolume() - pGeoData->innerVolume());
	// Compute isosurface normal vectors using function gradients
	calcSurfNormals(pGeoData);
	return pGeoData;
}

int CTreeContourExtractor::getStartTetra(const SuperArc& arc, float isoval)
{
	int v1, v2;		// order the vertices of arc such that val(v1) < val(v2)

	if (arc.intra_node < 0) // no intra node
	{
		if (m_pData->getValue(arc.v2) < m_pData->getValue(arc.v1))
		{
			v1 = arc.v2;
			v2 = arc.v1;
		} else {
			v1 = arc.v1;
			v2 = arc.v2;
		}


		if (!m_pVolCritical->areNeighbors(v1, v2)) {	// v1 and v2 are not neighboring vertices
			int xid;
			SuperNode& node1 = m_pConTree->getNode(v1);
			SuperNode& node2 = m_pConTree->getNode(v2);
			int x1 = node1.getNeighborTag(v2);
			int x2 = node2.getNeighborTag(v1);
			if (m_pData->getValue(x1) > m_pData->getValue(v1) &&
				m_pVolCritical->areNeighbors(v1, x1))
			{
				xid = x1;
				v2 = xid;
			} else if (m_pData->getValue(x2) < m_pData->getValue(v2) &&
				m_pVolCritical->areNeighbors(v2, x2))
			{
				xid = x2;
				v1 = xid;
			}
		}
	} else {
		v1 = arc.intra_node;
		v2 = m_pVolCritical->getCutEdge(isoval, v1);
		if (m_pData->getValue(v1) > m_pData->getValue(v2))	// swap v1 and v2
		{
			int tmp = v1;
			v1 = v2;
			v2 = tmp;
		}
	}
	int idx[3], nidx[3];
	m_pData->vert2index(v1, idx);
	m_pData->vert2index(v2, nidx);
	printf("idx : (%d %d %d) %f, nidx : (%d %d %d) %f\n", idx[0], idx[1], idx[2],
		m_pData->getValue(v1), nidx[0], nidx[1], nidx[2], m_pData->getValue(v2));
	int tid = m_pVolCritical->getTetraIndex(idx, nidx);

	return tid;
}

void CTreeContourExtractor::calcSurfNormals(GeometryData *pSurf)
{
	double grad[3];
	float *normals = new float[pSurf->numOfVertices()*3];
	for (int i = 0; i < pSurf->numOfVertices(); i++)
	{
		Vector3 v = pSurf->getVertex(i);
		m_pData->getVertGrad(v.x, v.y, v.z, grad);
		normals[3*i] = -grad[0];
		normals[3*i+1] = -grad[1];
		normals[3*i+2] = -grad[2];
	}
	pSurf->setNormalArray(pSurf->numOfVertices(), normals);
	delete[] normals;
}
